variables

ts=(0:stepsize:tend)'; %sampling times
N=length(ts);


testinp=@(t) sin(omega*t.*t); %input function
u=testinp(ts); %inputs at times


testt=tf([ucoeff;0]',[1;-ycoeff]',stepsize,'Variable','z');
% SYS = TF(NUM, DEN, TS, 'PropertyName1', 'PropertyValue1')
%creates a discrete-time transfer function with
%sample time TS (set TS=-1 if the sample time is undetermined), 
%with numerator(s) NUM and denominator(s) DEN.  
%The output SYS is a TF object.


ym=lsim(testt,u,ts);
%LSIM(SYS,U,T) plots the time response of the LTI model SYS to the
%    input signal described by U and T.  The time vector T consists of 
%    regularly spaced time samples and U is a matrix with as many columns 
%    as inputs and whose i-th row specifies the input value at time T(i).


y=ym;
firstentry=max(yl+1,ul);
ys=y(firstentry:N);


S = zeros(N-firstentry+1,length(coeff));



for k=1:yl
    S(:,k)=y(firstentry-k:N-k);
end

for k=0:(ul-1)
    S(:,k+1+yl)=u(firstentry-k:N-k);
end




disturbance=random('norm',0,ones(size(ys))*0.02);

yd=ys+disturbance;

plot(ts(firstentry:N),ys)




phat=S\yd;


